home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / minstd.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.8 KB  |  110 lines
  1. /* rng/minstd.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* MINSTD is Park and Miller's minimal standard generator (i.e. it's
  25.    not particularly good).
  26.  
  27.    The sequence is
  28.  
  29.    x_{n+1} = (a x_n) mod m
  30.  
  31.    with a = 16807 and m = 2^31 - 1 = 2147483647. The seed specifies
  32.    the initial value, x_1.  
  33.  
  34.    The theoretical value of x_{10001} is 1043618065, starting with a
  35.    seed of x_1 = 1.
  36.  
  37.    The period of this generator is 2^31.
  38.  
  39.    It is used as the RNUN subroutine in the IMSL Library and the RAND
  40.    function in MATLAB. The generator is sometimes known by the acronym
  41.    "GGL" (I'm not sure what that stands for).
  42.  
  43.    From: Park and Miller, "Random Number Generators: Good ones are
  44.    hard to find" Communications of the ACM, October 1988, Volume 31,
  45.    No 10, pages 1192-1201. */
  46.  
  47. static inline unsigned long int minstd_get (void *vstate);
  48. static double minstd_get_double (void *vstate);
  49. static void minstd_set (void *state, unsigned long int s);
  50.  
  51. static const long int m = 2147483647, a = 16807, q = 127773, r = 2836;
  52.  
  53. typedef struct
  54.   {
  55.     unsigned long int x;
  56.   }
  57. minstd_state_t;
  58.  
  59. static inline unsigned long int
  60. minstd_get (void *vstate)
  61. {
  62.   minstd_state_t *state = (minstd_state_t *) vstate;
  63.  
  64.   const unsigned long int x = state->x;
  65.  
  66.   const long int h = x / q;
  67.   const long int t = a * (x - h * q) - h * r;
  68.  
  69.   if (t < 0)
  70.     {
  71.       state->x = t + m;
  72.     }
  73.   else
  74.     {
  75.       state->x = t;
  76.     }
  77.  
  78.   return state->x;
  79. }
  80.  
  81. static double
  82. minstd_get_double (void *vstate)
  83. {
  84.   return minstd_get (vstate) / 2147483647.0;
  85. }
  86.  
  87. static void
  88. minstd_set (void *vstate, unsigned long int s)
  89. {
  90.   minstd_state_t *state = (minstd_state_t *) vstate;
  91.  
  92.   if (s == 0)
  93.     s = 1;    /* default seed is 1 */
  94.  
  95.   state->x = s;
  96.  
  97.   return;
  98. }
  99.  
  100. static const gsl_rng_type minstd_type =
  101. {"minstd",            /* name */
  102.  2147483646,            /* RAND_MAX */
  103.  1,                      /* RAND_MIN */
  104.  sizeof (minstd_state_t),
  105.  &minstd_set,
  106.  &minstd_get,
  107.  &minstd_get_double};
  108.  
  109. const gsl_rng_type *gsl_rng_minstd = &minstd_type;
  110.